
% Heaviside and its sensitivity
function [dH,dH_Phi]=NewHeaviside(phi,alpha,nnod,epsilon)
epsilon1=epsilon;
epsilon=1e-3;
num_all=[1:nnod]';
num1=find(phi>epsilon);
dH(num1)=1;
num2=find(phi<-epsilon);
dH(num2)=alpha;
num3=setdiff(num_all,[num1;num2]);
dH(num3)=3*(1-alpha)/4*(phi(num3)/epsilon-phi(num3).^3/(3*(epsilon)^3))+(1+alpha)/2;
dH=sparse(dH');

num1=find(phi>epsilon1);
dH_Phi(num1)=0;
num2=find(phi<-epsilon1);
dH_Phi(num2)=0;
num3=setdiff(num_all,[num1;num2]);
dH_Phi(num3)=3*(1-alpha)/(4*epsilon1)*(1-phi(num3).^2/((epsilon1)^2))+0.5;
dH_Phi=sparse(dH_Phi);dH_Phi=sparse(dH_Phi');
end